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Multivariate  Dependent  Renewal  Processes 


1.  Introduction.  The  standard  place  of  renewal  point-processes 
in  Reliability  theory  is  in  the  analysis  of  repair-vs-replace 
decisions  and  of  the  costs  associated  with  various  maintenance 
policies  (cf.  Barlow  and  Proschan,  1975,  Chapter  6).  However, 
even  when  attention  centers  on  a  single  device,  a  realistic 
model  of  reliability  would  require  description  of  dependent 
components  which  affect  each  other's  failure  rates  but  fail 
and  are  replaced  separately.  The  simplest  case  of  such  a  model 
would  be  a  system  of  independent  components  in  which  failure 


of  one  component  places  an  extra  load  or  shock  on  the  others, 
with  the  effect  of  permanently  multiplying  their  hazard  rates 
(by  a  fixed  function  of  time).  Freund  (1961)  introduced  this 
model  for  two  component  lifetimes  T^,T2,  defining  conditional 
hazard  functions  (cf.  the  definitions  in  Cox  and  Lewis,  1967) 

,  f  a  for  i  =  1 

lim  A  XP  { T .  5  t  +  A  |  min(T, ,T?)  >  t}  =  < 

A-*-0  [  p  for  i  =  2 


lim  A_1P{T.  5  t  +  A  |  T.  =  s  <  t  <  T .  } 
A-*0  1  3  1 


for 

i  =  1,  j  =  2 

for 

i  =  2,  1  =  1 

In  the  context  of  dependence  between  death  and  censoring 
times  in  survival  data  analysis,  Slud  and  Rubinstein  (1983) 
generalized  the  case  p  =  p'  of  the  Freund  model  to  allow  a, 
p,  and  a'  to  be  arbitrary  hazard  functions  depending  on  t. 

In  the  present  paper,  we  define  a  general  model  of  dependence 


2 


between  a  finite  number  k  of  component  lifetimes,  including  the 
Freund  and  generalized  Freund  models  as  special  cases.  The 
form  of  dependence  -  involving  the  current  life-lengths  of  the 
components  in  a  regression  model  like  Cox's  (1972)  for  the  con¬ 
ditional  hazards  -  allows  the  instantaneous  replacement  of 
failed  components  by  independent  new  ones,  leading  to  a  class 
of  dependent  failure  point-processes  generalizing  the  renewal 
point- processes . 

In  Section  2  we  define  our  model,  formulate  the  vector 
process  of  current  component  lifetimes  as  a  Markov  process  and 
study  the  behavior  of  stationary  distributions.  Some  explicit 
calculations  for  asymptotic  joint  distributions  of  component 
current- lives  are  given  in  Section  3  for  a  special  model  with 
one-way  dependence.  Also  for  two  generalizations  of  Freund's 
(1961)  shock  model,  we  discuss  in  Section  3  the  qualitative 
behavior  of  stationary  distributions.  Finally,  statistical 
inference  procedures  for  the  regression  coefficients  in  our 
model  are  derived  both  from  full  and  partial  likelihood  score 
statistics,  and  we  compute  asymptotic  relative  efficiency  of 
the  two  approaches  in  a  bivariate  example. 

2.  Model  assumptions.  Analysis  of  Markov  process.  We  suppose 
the  right-continuous  counting  processes  N^(t),  j  =l,...,k  are 
defined  on  t0,»)  and  have  unit  jumps  with  probability  1.  It 
is  natural  to  interpret  N^(t)  as  the  number  of  failures  (with 
instantaneous  replacement  by  an  independent  component)  up  to 
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time  t  of  the  jt^1  component  of  an  assembly.  Our  main  model 
assumption,  given  in  terms  of  conditional  hazard  functions  as 
defined  by  Cox  and  Lewis  (1972),  is 


(*) 


lim  h_1P(N.  (s+h)-N.  (s)  =  1  |  {N.(u):  0<uSs, 
h-0+  3  3  1 


=  (Tj (s) ) • exp( 


0ij<Pij(s(s)»  s 


1  <  i  <  k}) 
gj (x(s) ) 


where  for  s  >  0 


Tj(s)  =  s  -  min{u<s:  N^(u)  =Nj(s)} 
x(s)  =  (x^Cs) , . . . ,xk(s) ) , 


the  \ . ( * )  are  hazard  functions  and  lim  means  limit  in  the  mean. 
The  random  variable  x^(s)  is  the  current  life  at  time  s  of  com¬ 
ponent  j,  and  has  been  a  principal  focus  of  Reliability  theory. 
The  model  assumption  (*),  like  the  regression  models  of  Cox 
(1972)  which  inspired  it,  can  be  given  a  more  general  form  - 
the  log-linear  form  of  regression  is  especially  tractable  for 
statistical  inference,  but  by  no  means  necessary  for  the  Markov 
process  formulations  of  the  present  Section.  The  important 
feature  of  (*)  is  that  the  dependence  between  components  in¬ 
volves  only  their  current  ages  and  no  past  history.  (For  a 
model  where  past  cumulative  lifetimes  of  components,  from  the 
time  the  assembly  was  placed  in  service,  affect  current  hazards, 
see  Slud  and  Winnicki  (1982).)  We  remark  at  once  that  in  the 


simple  case  of  (ft)  where  all  are  constants,  it  is  clear 
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that  the  N  j  ( • )  are  k  independent  (delayed  unless  Tj(0)-=0) 
renewal  counting-processes. 

In  the  present  Section,  we  adopt  the  simpler  notations 


N 

gj(x>  =  xj(*j)exp<Z  Pij‘Pij(x))» 


lx)  =  f 
/  —  1 


gj<i) 


The  regularity  conditions  on  these  functions  needed  for  further 
work  are  contained  in  the  following  Theorem,  which  summarizes 
for  general  g j ( • )  the  information  available  from  standard  Mar¬ 
kov  process  theory.  Some  additional  notations  required  for  the 
Theorem  are: 


Ti(s)  =  inf  {t  >  s:  N(t)=N(s)+i} 


for  i  >  1 


A .. ( s )  =  j  for  which  Nj(T^(s))  -Nj(T^Cs)-)  =  1 

N(  • )  -  ^  N.  ( • ),  N(  ■ )  -  (Nj(  •),...  ,Nj^(  •)),  1_  -  (1, . . . ,  1)  €  ]R 

3*1  3 

=  (0,..., 0,1,0,..., 0)  =  the  j  canonical  basis  vector  in  B 

Xj  =  x  -  Xj  ej  for  x  €  1R 


Theorem  2.1.  Assume  that  are  right-continuous 

counting  processes  on  [0,«O  without  fixed  discontinuities, 
defined  on  some  probability  space  (&,f,P)  and  satisfying  (*) 
with  functions  gj ( • )  which  are  Borel  measurable,  non-negative, 
and  such  that  for  all  x 


fe  f 

j  g.(x  +  ul)du  i 


-*  “  as  e  -*• 


<  00  for  sufficiently  small  e=e(x)>0 
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for  j  =l,...,k,  where  e(x)  is  such  that  there  exist 
5  ’  eo  >  0  with  e(x)  >  eQ  whenever  min{Xj  :  li  j  <k}  <5. 

(i)  Then  fx(s):  s>0>  is  a  stationary-transition  Markov 
process  (with  piecewise-constant  right-continuous  paths) 
essentially  determined  by  the  joint  conditional  density  given 
Fg  =  o({N(u):  Oluis})  of  T-^(s)  and  A-^(s): 

f  u 

(2.1)  f ( u  ,  j  j  F  )  =  g.(x(s)  + ul)exp[-  I  g(x(s)  +vl)dv]. 

'o 

Moreover,  all  stationary  distributions  for  the  process 
{t(s):  s>0}  have  densities  n(-)  which  satisfy  for  Lebesgue 
almost  all  x 

k  r  ,  k  , 

(2.2)  Y.  J  n(Zj+xj£j)dZj  =  "  J  Tr(y)g(y;)dy  +  ]T  J  n(yOgj (y_)dy_ . 

[yi<xi,i?tj]  [yi<xi,l<i<k]  i  =  1  [yi<xi,i^j] 

(ii)  If  for  each  j  and  almost  all  x,  gj(x)  -  hj(Xj)  >0  for 

T 

measurable  h ^ ( •  )  with  /  g(x  +  ul)du  <  «  for  all  T < »  , 

«.  ru 

H.  r  I  exp(-H- (u))du  <  ”,  where  H.  (u)  =  I  h. (z)dz,  H.(+»)  =  »  ,  and 
3  ^0  3  1  -0  J  J 

g(x  +  ul)  -g(x!  +  ul)  >  hj  (Xj+u)  -  hj  (u)  for  x€[0,°°)\  u?0, 

then  the  laws  of  x(s)  converge  exponentially  fast  in  total 
variation  to  a  probability  law  with  density  n ( • )  uniquely 
characterized  by  (2.2). 
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Remark.  The  assumptions  (A)  turn  out  to  guarantee  that  the 
random  variables  T.j(s)  for  all  i  >  1,  s  >  0,  and  arbitrary 
t(0),  are  almost  surely  strictly  positive  and  finite. 

Proof  of  Theorem.  For  each  M  >  0 ,  s<u,  <u.<  ...  <uu<t, 
-  12  M 

M 

®tid  Ijkj  “  ^  {I>»**>k}  ,  we  define  u^  -  ( u^ ,  •  •  •  >  u^ ) 

and 

(2.3)  c^t,^,^  |  Fs)  =  lim  (2M61---6m)~1  •  P{N(s)  =N(u1-51), 

*»6M-*0+ 

N(ui+6i) -N(ui-6i)  =  e^  .,  l<i<M,  N(t)  = N(uM+6M) | Fg} 

where  existence  in  the  mean  of  the  limits  will  be  established 

inductively  in  what  follows.  Here  N( • )  is  almost  surely  right 

continuous  without  fixed  discontinuities  on  [s,°°)  given  F  ,  so 

s 

that  if  the  limit  in  (2.3)  exists  then  the  random  function 
CM^‘’— I  V  :’‘S  a,s*  continuous  on  [uM,«°).  In  particular, 
cQ(t|Fs)  =  P{N(t)  =N(s)  |  Fs}  exists  and  is  a.s.  continuous 
in  t  for  t?s.  Moreover,  if  we  denote  differentiation  from 
the  right  in  t  by  D+ ,  then  for  M>0  and  t>s,  existence  of 

(2.3)  together  with  (*)  implies 

(2.4)  D+cM(t,uM,iM  |  Fs)  =  -g(T(t))cM(t,uM,jM  |  Fg) 
where  T(t)  is  uniquely  determined  from  the  conditions 

Ti(s)=ui,  Af(  s )  =  j  i  for  1?  i?M;  TM+1  ( s )  >  t 


by  the  formulas 
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(2.5)  tCI^s))  =  T(Ti_1(s))  +  (Ti(s)  -  Ti_1(s))l 

-  (t.  (T.  ,(s))  +T.(s)  -T.  ,(s))e.  ,  2<iiH 

x(t)  =  x(Tr(s))  +  (t  -  Tm(s))1. 

From  (2.4)  it  follows  (Yosida,  1980,  p.  239)  that  a.s.  for 
t  >  s 

,t-s 

c0(t|Fs)  =  exp[-  J  g(T(s)  +  vl)dv] 

(2.6) 

rt-uM 

cM(t*1JM*WFs)  =  expC“  J0  S(1(V +vi)dv]cM(V  Fs) 


while  a  direct  application  of  (*)  to  (2.4)  implies,  assuming 
the  existence  of  (2.3)  for  given  M, 

(2.7)  CM+l^t’  =  ^  —  ^  ^CM  ^  ’^M  ^  Fs  ^  ' 


Now  (2.6)  and  (2.7)  show,  by  a  simple  induction  on  M,  that  all 
the  limits  (2.3)  exist  and  a.s.  for  M>0  (taking  uQ  =  s) 


(2.8)  y  ’  2m  ^  F  s  ^ 


m  u.-u.  . 

”  M  r  i  i-l 

gWu^J+vlJdv].  JTfexpC-j  g(l(ui_1)+vl)dv]g.  (t(u.))}. 

0  i— 10  1  i  ^ 


The  factored  form  of  the  functions  together  with  the  formu¬ 
las  (2.5)  determining  t ( • )  on  [s,t]  readily  implies  that 
{x(s)  ,  s  5 0 }  is  a  stationary-transition  (Strong-)  Markov 
process  on  [0,~)^. 


Now  we  need  to  make  some  technical  comments  regarding 


the  transition  functions 


P(t,x,A)  =  P{r(t )  (A  |  t(  0 )  =  x  } 
pn(t’*’A)  =  «  A  >  N ( t )  =  n  |  t(0)  =  x} 


where  A€B([0,*>)  )  is  a  Borel  set  and  n=(n^,n^ 
is  a  vector  of  nonnegative  integers.  First,  since 
tence  of  limits  in  (2.3)  implies  with  s = 0,  t(0) = 


the  exis- 
k 

M=  Y  n 
]  =  1 


(2.9) 


Pn(t,x,A)  = 


I 


/ 


CM(t,Uj, 
•<uM<t,  ; 


rA] 


(where  #  denotes  cardinality  of  a  set),  it  is  easy  to  check 
by  (2.8)  that  for  each  n,  x,  T 

X  /  /  pn<s,x,d^)g(y_)ds  =  P{N(T)  >  1  |  T(0)  =  x} 

n:n.>0, I5j5k  'o  J  2-  “ 

hence 


X  J  P  (t,x,dy^)g(^)  <  «  for  a.e.  t  (  [0,“). 

n:n^0,  lij^k  - 


Secondly,  from  formula  (2.9)  it  is  clear  that  if  > 1  for 
j=l,...,k,  then  pn(t,x, •)  is  absolutely  continuous  with  re- 
spect  to  Lebesgue  measure  on  [0,°°)  .  Next,  if  nj  -  2  for 
j  =  l,...,k,  (2.9)  implies 

(2,10)  ft  Pn(t’2i’CO’zl)x'*'xC0»V)  =  “  X  ferpn(t,x,C0,Z  )x...x[0, 

3  =  1  3  ~ 


f  k 

-  J  pn(t,x,dy_)g(£)  +  V 

CO,z1)x...x[o,Z]<;)  3  =  1 


/  Pn-£.(t»2i»^)gj(Z.) 
Cyi<zi,i^] 
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and  the  same  formula  holds  with  P  replaced  by  Z  Pn ■ 

-  n : n . >2  - 

Since 

J  g(x  +  u_l)du  =  «  for  all  x  implies 

yo 

Z  Pn(t,x,[0,»)k)  ►  0  as  t-» 

n:  nj<2  for  some  j  - 


the  last  statement  of  (i)  follows  easily. 

Finally,  in  part  (ii)  of  the  Theorem,  if  we  define  ran¬ 
dom  functions  cM  as  in  (2.8)  with  all  Cx>  replaced  by  hj(Xj), 
and  define  Qn  from  cM  just  as  Pn  was  defined  from  cM  in  (2.9), 
it  follows  inductively  from  assumptions  in  (ii),  for  fixed  M, 
0<u-L<u2<  ...  <uM<t,  and  j^,  if  r  ( *  )  is  given  on  [0,t] 
by  (2.5)  with  s  =  0,  that 


rt  c  t 

j  (g(x  +  vl)  -  g(m(v) ) )dv  5  I  (a(x  +  vl)  -  h(x(v) ) ) 


dv 


where  h(x)  =  V  h.(x.).  Therefore  (putting  s  =0  and  writing 

j^l  :  : 

cM( • >  for  CM( '  I Fn) 5 

exp  [J  g(x  +vl)dv]cM(t,uM,jM)  >  exp[  J  h(x  +  vl)dv]cM(t,uM,  j^) 

*0  0 


and  by  (2.9)  we  have 

■  t 


(2.11)  P  (t,x,A)exp[  f  g(x  +  ul)du]  tQ  (t,x,A)exp[  Y  (H .  (t  +  x. )  -  H.  (x . ) )  ] 
-  y0  -  j=l  3  J  J  1 


However,  Qn(t,x,A)  are  precisely  the  transition  functions 
P{t^( t )  (  A  ,  N^( t )  =  n  |  t( 0 )  =  x  }  for  a  process  ( ( s ) , ( s ) :  s  t  0 ) 
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satisfying  (*)  with  (x)  replaced  by  hj(Xj),  that  is,  for  the 
current-life  and  counting  processes  associated  with  k  inde¬ 
pendent  renewal  processes  with  renewal  hazards  h ^ ( • ) .  The 
assumption  M-j  <  ®  implies  these  renewal  current-life  processes 
are  positive-recurrent  with  (Karlin  and  Taylor,  1975,  p.  193) 


(2.12)  Q ( t ,x, A) 


X  Qn<t,x,A) 
n>0  - 


exp[- 


t 


H j (y j ) ]dy 


Q(A) 


as  t  -*■ 

(2.13) 


■> ,  uniformly  in  x  €  [0,“)  ,  A  f  B(tO,»)  )•  But 

k  r  t 

P(t,x,AC)  >  Q(t,x,Ac)exp[  y  (H.(t+x.)  -H.(x. ))  -  I  g(x+ul)du] 
-  -  ^71  J  J  J  J  0 


implies,  since  P(t,x, •)  and  Q(t,x, •)  are  probability  measures. 


(2  .14)  P( t ,x, A)  < 
k 


1-expCt"  (H.(t+x.  )  -H.(x.))  -  f  g(x+ul)du]-(l-Q(t,x,A)). 
JFl  3  3  3  3  y0 


Fix  1  >  5  >  0 .  By  comparison  with  the  independent  renewal 
processes  with  hazards  h^ ,  one  can  find  tg  so  large  that 

(2.15)  P  {Nj  ( tQ)  >1  for  j  =  1,  .  .  .  ,k  |  t(  0 )  =  x  }  1-6 

uniformly  in  x€[0,»)  .  Then  (2.12)  says  that  for  ^  large 
enough,  for  all  x  and  A, 


(2.16) 


|Q(t1-t0,x,A)  -  Q( A) |  <  6. 


1 

■! 


4 
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On  the  event  CNj(t0)>l  for  j  =  l,...,k],  by  local  integrabil- 
ity  of  g,  there  exis*-  61>0  such  that,  for  the  fixed  tx 
in  (2.16) 


(2.17)  inf 

t(0)=x*[0 


* 

,.)  ^j'vvyv' -Hj(''j<to>)) 


-P 


g(x(t0)  +ul)du]  >  6.  . 


Therefore  by  (2.14)  -  (2.17),  for  all  x  6 [0,~)K  and  all  A 

(2.18)  P(t1,x,A)  =  J P(tQ ,x,d^)P(t1-t0,^,A) 

<  6  +  (1-S)[1  -  61  +  6,(6  +  Q( A) ) ] 

=  1  -  6. (l-5)tl  -  6  -  Q( A) ] 


so  that,  with  e  s  min((l-6)/2  ,  6^1-6)  /2)) 

(2.19)  if  Q(A)  i  e,  then  P(tl5x,A)  <  1-e. 

This  is  Doeblin's  condition  (D)  (Doob,  1953,  p.  256)  and  by  a 
Theorem  proved  in  Doob  (1953)  our  Theorem  is  established. 

By  the  method  of  characteristics  (Courant,  Hilbert,  1962, 
pp.  170  ff.)  it  follows  immediately  from  (2.2)  in  Theorem  2.1  (i) 
that  for  x  €  [0,°°)k 


(2.20) 


f  *3 

rr  (x)  =  rr  (x-x.l)-exp[-l  g(x  +  (u-x- )l)du] 

"  3“  /O  J 

ir(x-Xjl)=  j*  tt  (x-Xji+ye^  )gj  (x-x  -l+ye.  )dy 


when  x.  =min  x-. 
■*  i 
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The  system  (2.20)  of  integral  equations  is  already  a  non¬ 
trivial  generalization  of  the  renewal  equations  which  determine 
the  current-life  asymptotic  density  «(•)  in  case  all  0^  =0 
(see  Karlin  and  Taylor,  1975,  pp.  192-3).  Theorem  2.1  contains 
useful  information  on  existence  and  uniqueness  of  solutions  of 
systems  of  the  general  form  (2.20). 

To  apply  these  asymptotic  results  in  reliability  theory, 
one  would  want  also  to  know  the  asymptotic  distributions  of  the 

times  between  successive  jumps  for  each  of  the  counting  processes 

k 

N. ( • )  and  for  the  superposition  counting-process  )  N.(')SN(>). 

J  jTl  3 

These  are  easy  to  derive  by  ergodicity  of  the  process  {t(s)> 
with  initial  density  n(-),  using  a  well-known  generalization 
(e.g.,  see  Rolski,  1981,  Ch.  1)  of  a  Strong  Law  technique  given 
by  Karlin  and  Taylor  (1975).  We  state  only  the  result  for  den¬ 
sities  of  single  times  between  replacements,  although  corre¬ 
sponding  asymptotic  joint  densities  of  successive  times  between 
replacements  can  also  be  readily  expressed  in  terms  of  tr(  •  ) . 

Proposition  2.2.  Under  the  hypotheses  of  Theorem  2.1  ensuring 
the  validity  of  (2.2),  if  (2.2)  has  a  unique  solution  and  we 
define 

W^m)  s  min{ t  ^  0 :  N  ^  ( t )  =  m  }  -  mini t  >  0 :  N.  ( t )  =  m-1 } 

W(m)  s  min{  t  >0:  N(t)  =m>-  min  { t  >  0 :  N(t)  =  m~l } 

then  as  m  ■* «  the  asymptotic  densities  f  ( •  )  of  and  f  ^  (  • ) 

of  W<m)  on  [0,”)  are 
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(2.21) 


vx)  =  -a*  /  ><«*;//  "<^>^3 

[o,-)*"1  [o,-)k  1 

f(x>  .  -fe  /  |  nt^)^// 

[x,~)  J  LU,®> 


Corollary  2.3.  The  asymptotic  marginal  densities  (under  ( ii) 
of  Theorem  2.1)  as  s  -» ■» 

J  n (y’! +xe . )dy/.  for  x.(s) 

'  T,  -I  3  3  3  J 

r  n  \  K  -L 


[0,0°) 


and 


)  J  rr(y!+xe.)dy!  for  min  {t,  ( s ) , .  .  .  ,x.  (s )  } 

t=s  '  i,_i  3  -]  -3 

3-1  [0,»)k  1 

are  non-increasing  functions  of  x. 

3.  Examples.  Consider  the  model  of  the  previous  Section  in 
the  case  k  =  2,  <p21(*)»0,  P11  =  P22~0,  P12  =  1’  X2('^*X2 

constant,  and  <p12(x,y)=ln  h(x).  Then  g1(x,y)  =  X^(x) , 

g  (x,y)  =  h(x)X  ,  and  we  define  H(x)  =  f  h(z)dz,  A,(x)  = 

0  l 

fX 

J  X-^(z)dz.  Since  g-^(')  depends  on  x  alone,  the  counting 

process  N^t*)  will  marginally  be  an  ordinary  (delayed)  re¬ 
newal  counting  process  with  interoccurrence  distribution  func- 

fX 

tion  F.(x)  =  1  -  exp(-  )  X.(z)dz),  while  N0(*)  will  be  a 

'  0  i 

"doubly  stochastic  Poisson"  process.  We  will  calculate  for 
this  example  an  explicit  form  for  the  asymptotic  joint  density 
rr  ( •  )  of  current  lives  t^(s),t2(s). 
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First,  we  observe  that  the  renewal  counting  process  N^(*) 
has  corresponding  stationary  distribution  if  and  only  if 

(3.1)  p-,  =  /  exp(-A,  (t )  )dt  <  ». 

1  'o  1 

Then  the  condition  P{N(t)<-}  =  1  for  all  t  is  equivalent  to 

P{N^(t)=0,N2(t)<«}  =  exp(-A^(t))  for  all  t  >  0 

which  says  that  with  probability  1  there  are  finitely  many  jumps 
in  N^(*)  between  successive  zeroes  of  t^(*).  But  between  zeroes 
of  ^(O  is  the  counting  process  for  a  nonhomogeneous 

Poisson  process  with  rate  X^hC*).  Therefore 

(3.2)  P{N(t)<«)  =  1  for  all  t  if  and  only  if 
H(t)<~  for  all  finite  t  with  A^(t)<«. 

Under  assumptions  (3.1)  and  (3.2),  it  is  easy  to  check  that  the 
hypothesis  of  Theorem  2.1  (ii)  is  equivalent  to  as  ess.inf.h(’) 

>  0,  together  with  A^(“)  =  A-^(x)  +  H(x)  <  ~  when  x  <  <». 

Rather  than  impose  these  last  hypotheses,  we  calculate  directly 
the  unique  solution  n(x,y)  of  (2.20).  Theorem  2.1  (i)  then 
implies,  by  standard  arguments,  that  {t(s):  s>0)  is  ergodic 
and  for  all  initial  distributions,  as  s  -*■  -  the  distributions 

2 

of  x ( s )  converge  weakly  to  the  measure  with  density  a  on  [0,~)  . 


Equations  (2.20)  can  be  rewritten  in  our  example  to  give 


(3.3)  n(x,y)  = 


n(x-y  ,  0)expC-A1(x)  +  A1(x-y)  -  X2(H(x) -H(x-y))]  if  x>y 


”(0,y-x)exp[-A1(x)  -  X^H(x)] 


if  x<y 


n(x,0)  =  /  X  h(x)rr(x,y)dy,  n(0,y)  =  /  X,(x)n(x,y)dx. 

/  r*  *-  •  n 


Since  N^( • )  is  a  renewal  counting  process  with  interoccurrence 
time  distribution  function  1  -  exp(-A ^x) ) ,  standard  renewal 


theory  implies 


and  by  (3.3), 


/  "(x. 


y)dy  =  p~1exp(-A1(x) ) 


n(x,0)  =  p11X2h(x)exp(-A^(x) ) 


(3.4) 


n(0,y)  =  J  X1(x)n(0,y-x)expC-A1(x)  -  X2fi(x)  ]dx 
C  “ 

+  j  X2n.^1X1(x)h(x-y)e:<p[-A1(x)  -  X2(H(x)  -  H(x-y))]dx. 


From  the  second  of  these  equations  follows 


oo 

(3.5)  q(t)  =  J  exp(-ty)n(0,y)dy 
0 

rm  X,(x)  -A,  (x)-X0H(x)-t>:  r x 

=  J  -i -  e  1  2  r  1 

\  “l 


/■Jf  X  H(z)+tz 

[  J  X„h(z)e  1  ' 

0 


dz]dx 


/«  -rx-A, 

-  A  vx>e 


-tx-A,(x)-X0H(x) 


Then  the  known  Laplace  transform  q(t)  determines  it(0,y),  and 
(3.3),  (3.4)  give  the  exact  asymptotic  joint  density  n(x,y). 
In  particular,  the  form  of  (3.3)  immediately  implies  that 


\ 
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under  the  model  of  this  Section  as  s-*-«,  for  any  initial 
distribution  irQ(x,y), 

t1(s)  and  i:2(s)  -  t1(s)  are  asymptotically  conditionally 
independent  given  sgn<T2 ( s )  -  t^( s ) ) . 

There  are  two  generalizations  of  Freund's  (1961)  bivariate 
shock  model  which  may  be  of  interest  in  Reliability.  Let  k  =  2, 

=  ?22  =  anc*  =  ^21  =  ^  ^  ^  *  Then  the  two  general 

models  are  given  by 

<P12(x,y)  =  Ij-  3  In  hx(y) 

(pi> 

[  «)21(x’y>  =  I[ycx]  ln  h2(x) 

and 

<P12(x,y)  =  If  3  In  h1(y-x) 

(f2) 

_  (p2i(x»y>  =  I[y<x]  ln  h2(x"y) 
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where  1^  j  denotes  indicator  function.  In  these  models,  com¬ 
ponent  1  acquires  the  multiplicative  factor  h2(*)  in  its 
hazard  the  instant  that  component  2  fails;  however,  further 
failures  of  the  replacements  of  component  2  during  the  life 
of  the  same  component  1  do  not  further  change  the  hazard. 
Model  ( F-^)  says  essentially  that  failure  of  one  component  in¬ 
stantaneously  shifts  the  other  component  to  a  new  survival 
curve,  while  (F?)  says  that  the  proportionality  factor  for 
hazard  due  to  shock  can  vary  with  time  only  as  measured  from 
the  instant  of  shock.  The  bivariate  model  of  Slud  and  Rubin¬ 
stein  (1983)  is  a  special  case  of  (F^). 

Under  model  (Fx),  (2.20)  easily  implies  the  asymptotic 
density  for  x  has  the  form 


j  n(0,y-x)exp[-A1(x)  - 
(3.6)  n(x,y)  =  < 

(y 

tt  (x-y  ,0)exp[-  J  h2( 


so  that  under  (F^) 


j  h  (y-x+u)X  (y-x+u)du] 
0  1  l 


x-y+u)X1(x-y+u)du  -A2(y)] 


if  x  < 
if  x  > 


(3.7)  min(x1(s),T2(s))  and  max(x^(s),x2(s)) -minute), x2(s)) 

are  asymptotically  conditionally  independent  given 
sgn(x^(s) -x2(s))  as  s  -*■». 


Although  no  such  appealing  property  can  hold  for  nr  ( •  )  under 
(F?),  except  in  Freund's  case  with  X.  and  constant,  it  seems 
to  be  physically  the  more  reasonable  of  the  two  models  and 
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should,  together  with  the  statistical  procedures  of  the  follow¬ 
ing  Section,  find  applications  in  Reliability.  It  is  interest¬ 
ing  also  to  remark  that  if  in  (F^)  one  has  h^(*)  sharply  peaked 
at  0  wit'n  values  rapidly  decaying  to  1,  then  model  (F2)  mimics 
the  behavior  of  the  well-known  model  of  Marshall  and  Olkin  for 
dependent  failure  times  (see  Barlow  and  Proschan,  1981,  Chap¬ 
ter  5). 


4.  Statistical  inference  for  dependent  renewal  processes. 
Until  now,  we  have  made  no  use  of  the  special  regression  form 
for  g j ( •  )  in  (*).  We  recall 

g.(x)  =  X  .  ( x .  )exp(  B..<j>..(x)) 

J  ~  J  J  iJ  iJ  ~ 


which  for  convenience  (and  identif iability  of  parameters  p^.) 
we  rewrite 

g,.(x)  =  \j(x)exp(£.  PijVijCx)). 

Now  suppose  that  data  on  the  process  (x(s)  ,  0  <s  5t}  is  avail¬ 
able  in  the  following  form: 


t( 0 )  and  {T^,A^:  £=l,...,N(t)}  are  observed 

where 

Tt  =  inf { s  >  0  :  N(s)  =  6},  =  j  iff  t^(T^)=0. 

The  full  likelihood  L  =  L(t(  0 ) ,  (T^  ,  ) ;  (p„  ) )  conditional  on 

t(0)  for  this  data  is 

N(t)  k  !U,  =  j]  At 

L  =  IT  "IT  Cg -  <'r<T-~)>3  •  exp(  -  J  g(t(u))du) 

1*1  j=l  3  ~  1  0  _ 
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which  makes  sense  because  for  all  l  ?0  (with  Tg  s  0 ) 


I'W  *  l(Te>  *  <Te  +  i  "  Te)i' 


The  score  statistic  for  P^.  is  given  by 

<4-1)  WTT1061-  =  ^  -  (  *y(l<u»gj(l(u))du 

from  which  the  identif iability  of  p^  for  each  j  and  i  t  j  is 
seen  to  be  assured  if  the  known  measurable  functions  <j>^  for 

i  i  {1, _ ,k)\{j  }  are  linearly  independent  on  a  subset  with 

positive  Lebesgue  measure  of  {x:  \Hx)  >0}.  The  log-concavity 
of  L  follows  from 

a2  (X 

(4.2)  -5 — — - log  L  =  -Ip.f_.-J  I  •  •  (t(u) )<p .  i  ■  (r(u) )g. (t(u) )du. 

'x  j'  L3  -3J  13  ~  i  3  ~  3  ~ 


If  we  assume  that  (2.2)  has  a  unique  solution  n,  so  that 
(t(s),  s  >0}  is  an  ergodic  Markov  process  by  Theorem  2. l(i),  then 
for  (P-^jZ  i,j=l,...,k,  i  *  j  )  in  an  open  region  D  c  j^kCk-l)  ^ 


as  t  -*■  ■ 


(4-3)  ~  h— a«-T  ;  log  L  ~  ti  ,  /  <p..( 


ap.-ap.,., 


/  (Pi^(x)<Pi'j(x)gHx)”(x)dx 

[0,-)k 


Moreover,  for  large  t  the  maximum- likelihood  estimates  (5--) 

^13 

for  which  =  0  are  unique,  strongly  consistent,  and  asymp¬ 

totically  jointly  normal  with  asymptotic  variance-covariance 

p  p 

matrix  the  inverse  of  the  matrix  tip  with  ij,i'j'  component 
equal  to  the  right-hand  side  of  (4.3)  (see  for  example  Basawa 
and  Prakasa  Rao,  1980). 
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Another  framework  for  the  consistent  estimation  of  (B..) 
is  given  by  the  recent  generalization  by  Andersen  and  Gill 
(1982)  of  the  partial  likelihood  methods  of  inference  due  to 
Cox  (1972,  1975).  Following  Cox  and  Andersen  and  Gill,  we  de¬ 
fine  the  partial  likelihood  for  ( )  based  on  the  data 
{T l  =  1, . . . ,N(t)  } 

(P)  N(t )  /  k  CA»  =  j]\  / 

(4.4)  =  TT^  (  TT^  gj(i(Tr))  j/g(l(Tr)) 


leading  to  the  "partial- like lihood  score  statistic  for  " 


(4.5) 


3P 


ID 


(P)  f  g^TCr.-))  1 


all  linear  combinations  of  which  are  martingales  in  the  con¬ 
tinuous  parameter  t.  As  in  the  maximum-likelihood  theory,  an 
argument  exactly  analogous  to  that  in  Andersen  and  Gill  (19 82) 
shows  that  for  large  t  the  maximum  partial  likelihood  estimates 

(8..  „)  where  (4.5)  equals  0  are  unique,  strongly  consistent, 
KiD  ’ P 

and  asymptotically  jointly  normally  distributed  with  asymptotic 
variance  the  inverse  of  the  matrix  with  ij?i  D  entry 


(4.6)  - 


0pijapi,Df 


log  L 


(P)  _ 


N(t) 


I[j  =  3 ' 3 


g.(x(T-)) 

(pij(-(V))<pi'j(-<V))  g(T(Tr)) 


g.(r(T  -)) 
1  “  g?T(Tr)l 


The  right-hand  side  of  (4.6)  is  asymptotic  in  probability  as 
t  -+ «  to  a  quantity  til?)  •  •  •  »  •  »  which  can  be  explicitly  calcu- 

lated  in  terms  of  n  in  special  cases. 
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We  mention  the  partial  likelihood  methods  of  estimation 
so  prominently  for  several  reasons.  First,  with  their  most 
recent  j ustif ication  by  Andersen  and  Gill  (1982),  they  provide 
easily  constructed  consistent  estimates  in  many  cases  where 
ordinary  maximum  likelihood  methods  are  not  available:  for 
example,  where  all  gj (x)  contain  a  common  but  unknown  "nui¬ 
sance"  factor  \(x),  or  where  data-collection  is  subject  to 
independent  random  censoring  with  unknown  distribution.  (We 
give  no  detailed  discussion  of  either  of  these  situations.) 
Another  reason  for  introducing  partial  likelihood  methods 
here  is  the  broad  exposure  they  have  been  given  in  the  recent 
literature  on  statistical  analysis  of  failure  time  data  (see 
for  example  Kalbfleisch  and  Prentice,  1980,  Miller,  1981). 

That  literature  showed  for  independent  failure  times  the 
considerable  theoretical  and  practical  efficiency  of  partial 
likelihood  methods  relative  to  maximum  likelihood.  In  the 
present  context  of  dependent  failure  times,  we  show  by  calcu¬ 
lations  along  the  line  of  Slud  (1982)  in  an  example  that  par¬ 
tial  likelihood  estimation  may  typically  be  quite  inefficient 
in  ignoring  the  lengths  of  intervals  between  events. 

Consider  the  model  (*)  with  k=2,  \^(x)  general, 

\?(x)  =  \2  constant,  P11  =  P22  =  P21  =  0,  p12  =  p,  and  q>12(x,y) 

Ir  t  Again  the  condition 
Lx<yJ  & 

=  J  exp(-A.  (x))dx  <  °° 

1  y0  1 
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guarantees  that  N^(*)  is  the  counting  process  for  positive- 
recurrent  renewals.  By  (2.5)  we  find  the  unique  solution  of 
(2.3)  satisfies 


"(x,y;P)  =  n  (  0  ,y-x;P  )exp[  -  f  ( X  ( u )  + ePX  >du ]  for  x<y. 
By  the  known  asymptotic  marginal  behavior  of  t-^(s), 

/oo 

n(x,z‘,p)dz  =  p...  exp(-A . (x) ) ,  x>0. 

N  XX 


Therefore  the  quantity  7*  from  (4.3)  is  given  in  the  present 
example  by 

(4.7) 


ll  ~  ^  e^X9rr(x,y;p)dx  dy 

P  Cx<y ]  2 


Lx< y  ] 


e^X^Tt  ( 0  ,y-x ; p  )  •  exp(-A ^(x)  -  e^X^x)dy  dx 


=  ePX^  J  expC-A1(x)  -  e^X^xldx  j  Tr(0,z;p)dz 


0 

p.  -i  r 

*  Vl  A 


0 


exp[-A1(x)  -  epX^x]dx. 


From  (4.6)  we  find  that  h  is  given  by 


I  =  p.lim  t 


-1^,  .  ...  VT2(Vm2e 


T^V)<T2(Te")3 


.M 


(X  (T2(Te-))  +  V  ) 


Now  a  failure  with  l  >2  can  have  t-^T^-)  < 

only  if  it  follows  a  failure  with  A<(_1  =  1.  Moreover,  immedi¬ 
ately  following  the  failure  at  T^_1  with  A^_1  =  1  the  con¬ 
ditional  hazard  for  N2(-)  is  the  constant  X2exp(p),  so  that 
the  time  from  until  the  next  jump  in  N2(-)  is  exponential 
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with  parameter  X^expCp),  and  the  time  has  condi¬ 

tional  survival  function  expC-X^e^  -A^(x)).  Since  the  number 
N1(t)  of  failures  with  =  1  is  asymptotic  to  t/^  by  the 

elementary  renewal  theorem,  we  conclude 


(4.8) 


=  vV  /, 


P 


o 


Vy)  e 

- o  exp ( -X  er  -  A, (y  )  )dy . 

X1(y)+X2eP  1  X 


p  r 

From  (4.7)  and  (4.8),  it  is  apparent  that  7^  <  for  all 

p,  X2,  X2(.),  with  the  interpretation  for  X^(*)  uniformly 
bounded  by  a  constant  X^  that  the  asymptotic  relative  effi¬ 
ciency  of  maximum  partial  likelihood  to  maximum  likelihood 

estimation  of  p  in  this  example  is  at  best  X^7  ( X"  +  X2e^), 

P  F 

and  I  „  / 1  a  can  be  nearly  1  only  when  X  (■)  is  very  large, 
p  p  1 

The  special  forms  of  the  estimators  may  also  be  of  interest 
here : 

N<i)  ,«Lt) 

A  (I  ^ue-2 , 


•  \  t  A 

where  Ag  =  1  iff  t-^(O)  <  t2(0),  while  Pp  is  the  unique  solu¬ 
tion  p  of 


0 


n  ( t )  r 

i  (  I[  V2  ’  A«-l=1] 


V 


[A«-1=1]  X2ep  +  X1(t1(T^-)) 
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